Momentum distribution of a freely expanding Lieb-Liniger gas 
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We numerically study free expansion of a few Lieb-Liniger bosons, which are initially in the 
ground state of an infinitely deep hard-wall trap. Numerical calculation is carried out by employing 
a standard Fourier transform, as follows from the Fermi-Bose transformation for a time-dependent 
Lieb-Liniger gas. We study the evolution of the momentum distribution, the real-space single- 
particle density, and the occupancies of natural orbitals. Our numerical calculation allows us to 
explore the behavior of these observables in the transient regime of the expansion, where they are 
non-trivially affected by the particle interactions. We derive analytically (by using the stationary 
phase approximation) the formula which connects the asymptotic shape of the momentum distribu- 
tion and the initial state. For sufficiently large times the momentum distribution coincides (up to a 
simple scaling transformation) with the shape of the real-space single-particle density (the expansion 
is asymptotically ballistic). Our analytical and numerical results are in good agreement. 

PACS numbers: 05.30.-d, 03.75.Kk, 67.85.Dc 
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I. INTRODUCTION 

Exactly solvable models describing interacting bosons 
in one-dimension (ID) have been studied over decades 
since thepioneering work of Girardeau [l[ , and Lieb and 
Liniger [3| . The interest in these models is greatly stimu- 
lated with recent experiments [1, 0, @ , in which ultracold 
atomic gases arc tightly confined in ID atomic waveg- 
uides, such that transverse excitations are suppressed. 
The Lieb-Liniger (LL) model describes ID bosons with 
pointlike contact interactions of a given strength c 
In the limit of sufRciently strong interactions, the LL 
gas enters the Tonks- Girardeau (TG) regime of impen- 
etrable bosons the TG regime can be obtained at 
very low temperatures, with strong effective interactions, 
and low linear particle densities p, 0, [l] • An interesting 
aspect of ID Bose gases, which can be probed experi- 
mentally from weakly to the strongly interacting regime, 
is their behavior out of equilibrium (e.g., see Ref. [H])- 
An exact (analytical or numerical) theoretical calculation 
of nonequilibrium dynamics of a LL gas is a complex 
many-bod y p roblem, which was studied in a few cases 
[E) m, 111. Il2l. A paradigmatic problem in this con- 
text is one-dimensional free expansion from an initially 
localized state [H [il [H, H [13, [3 . Quite generally, 
free expansion can be used to provide information on the 
initial state (e.g., see Refs. [1^ [13, [Sj and references 
therein) . 



Free expansion of a LL gas was studied in Ref. [14 1 
by using the hydrodynamic approximation it was 
demonstrated that the evolution of single-particle (SP) 
density is not self-similar for finite interaction strength 
c. Exact solutions of free expansion were studied both 
in the TG regime (c = oo) [m[llE3,[ll, and for a re- 
pulsive LL gas (0 < c < oo) [ij, [l3|. In the TG regime, 
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it has been shown that the momentum distribution ap- 
proaches that of noninteracting fermions during free ex- 
pansion [Tsl . [l6t . This was shown numerically by using 
the model of hard-core bosons on the lattice in Ref. [31 , 
and by using the stationary phase approximation in the 
continuous TG model (for the initial harmonic confine- 
ment) in Ref. [l^ . For quite general initial conditions, 
the asymptotic form of the wave functions for a freely 
expanding LL gas was calculated in Ref. (isj by using 
the stationary phase approximation (free expansion for 
a particular family of initial conditions was previously 
considered in [l3l)- R was shown that these wave func- 
tions vanish at the hyperplanes of contact between par- 
ticles, which is characteristic for TG wave functions [l|. 
However, it was emphasized that the properties of such 
asymptotic states can considerably differ from the phys- 
ical properties of a Tonks- Girardeau gas in the ground 
state of some trapping potential [l3| (see also the sec- 
ond item of Ref. [i3l)- We also point out that expan- 
sion dynamics is subject of studies in other ID models; 
for example particles interacting via an inverse-square 
pair potential , strongly correlated fermions [1^ , and 
strongly interacting Bose- Fermi mixtures [IJ]. To the 
best of our knowledge, studies of the momentum distri- 
bution of an expanding LL gas, which are based on exact 
time-dependent solutions, have not yet been made. 

Here we numerically study free expansion of a few Lieb- 
Liniger bosons, which are initially in the ground state 
of an infinitely deep hard- wall trap. The numerical cal- 
culation is carried out by employing a standard Fourier 
transform, as follows from the Fermi-Bose transformation 
for a time-dependent Lieb-Liniger gas d, [l3| . We focus 
on dynamics of one-body observables of the system, in 
particular the momentum distribution, the occupancies 
of natural orbitals, and also the real-space single-particle 
density. Our numerical calculation allows us to explore 
the behavior of these observables in the transient regime 
of the expansion, where they are non-trivially affected by 
the particle interactions. We derive analytically (by using 
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the stationary phase approximation) the formula which 
connects the asymptotic shape of the momentum distri- 
bution and the initial state. For sufRciently large times 
the momentum distribution coincides (up to a simple 
scaling transformation) with the shape of the real-space 
single-particle density, reflecting the fact that the expan- 
sion is asymptotically ballistic. The relation between the 
asymptotic expansion velocity of the LL cloud, and the 
overall energy stored in the system is derived. Our ana- 
lytical and numerical results are in good agreement. 

Before proceeding, let us devote a few words to some 
of the techniques for solving LL and TG models. The 
Bethe ansatz can be used to find the eigenstates for LL 
particles on an infinite line 0], with periodic boundary 
conditions 01, and in an infinitely deep box [2^ . In the 
TG limit; both stationary [l| and time-dependent [26j 
wave functions, in an arbitrary external potential, are 
constructed by using the Fermi-Bose mapping, i.e., by 
solving the Schrodinger equation for spinless noninter- 
acting fermions, after which the fermionic wave func- 
tion is properly symmetrized to describe TG bosons 
P, [2^. In a similar fashion, exact time-dependent 
LL wave functions, in the absence of an external po- 
tential and on an infinite line, can be constructe d by 
employing the Fcrmi-Bose mapping operator Il2l |. 
which ensures that the so-called cusp-condition (see e.g., 
[lOl . IT^ ) imposed by the interactions is obeyed during 
time-evolution. In the light of the recent experiments 
d, H, Q , there has been renewed interest in exact studies 
of LL gases ^stationary [13, [H, [H, HO, [Ml, [H and time- 

and time-dependent TG gases 



dependent Il0l. ll 

Even once a wave function describing a LL gas is 
known, determination of its correlation functions and ob- 
servables such as the momentum distribution is difficult 
m. Hi, EO, Ell, m, m. El, m. various methods were 
developed over the years including the Quantum Inverse 
Scattering Method [H, EH| j 1/c expansions [3§| from the 
TG (c 00) regime, and Quantum Monte Carlo Integra- 
tion [4^ . In the TG limit, momentum distribution can 
be analytically studied in a few cases (e.g., see 0, E3)- 
Numerical calculation of the TG momentum distribution 
can be performed also for excited and time-dependent 
states [3, [131; for hard-core bosons on the lattice see 
[l5| ; a simple formula suitable for numerical calculations 
was recently derived for the continuous TG model [s^ ]. 
and generalized for ID hard-core anyons [48|. 



II. THE LIEB-LINIGER MODEL AND 
OBSERVABLES OF INTEREST 



We examine a system of identical J-interacting 
bosons which are constrained to one spatial dimension. 
The Schrodinger equation for that system, in the absence 



of any external potential, is 



_ O ipB 

dt ~ ^ dx^ 



J2 2cS{x,~Xj)^B- (1) 

l<i<j<N 



The strength of interaction is described by a parameter 
c (here we consider repulsive interactions c > 0). The 
initial condition 'ijjBixi,X2, ■ ■ ■ ,XN,t = 0) is a localized 
state, e.g., the LL ground state in some external trapping 
potential. Since we are dealing with symmetric (bosonic) 
wave functions, it is convenient to write Eq. ([T]) in one 
permutation sector of the configuration space, Ri : xi < 
X2 < ■ ■ ■ < Xn, 



dt ^ dxf 



(2) 



the (5-function interactions are equivalent to the bound- 
ary condition [3|, 



1 - 



d 



dx 



_d_ 

dxi 



V'B =0. 



(3) 



It is worth to mention that the so-called cusp condition 
in the TG limit (when c — > 00) reduces to the condition 
that the wave function vanishes whenever any of the two 
particles touch. 

The time-dependent Schrodinger equation ([T]) can be 
solved exactly by employing the Fermi-Bose mapping op- 
erator [i,[ii|. 



n 

l<i<j<N 



sgn(Xj 



1 
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\dxj 


dxi J 



(4) 



If we find a fully antisymmetric (fermionic) wave function 
ipF which obeys 



.dipt 



E 



dt ^ dx'j 
then the wave function 



(5) 



(6) 



where A/'c is a normalization constant, obeys Eq. ([T]) 

i[i3,[il. 

This means that free expansion solutions can be found 
by using Fourier transform. Let 



iPf{xi, . . .,XN,t) 



dki ■ ■ ■ dk 



N 



where uj{k) = fc^, denote an antisymmetric wave function 
TpF, which evidently obeys Eq. ([5]); here 



ipFiki, ■ ■ ■,kN) 



1 



dxi ■ ■ ■ dx 



(27r)^ 

X VF(a;i,...,a;Ar,t = 0)e"'^j"=i''- 



N 



(8) 
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is the Fourier transform of ipp a,t t = . By acting 
with the Fermi-Bose mapping operator on ipp [see Eq. 
([6|)] wc obtain a time-dependent wave function of a freely 
expanding LL gas: 



ipB,c{^i, • ■ ■ , XN,t) ^ J dki-- ■ dkN 



(9) 



where 



G(A:i, . . . ,A:jv) = Nc^piki, ■ ■ ■ ,kM) 



X n 

l<i<3<N 



s.gn{xj - Xi) + - {kj - ki) 



(10) 



The information on initial conditions is contained in 
G(fci, . . . , fc^v) [that is, within . . . , /cAr)]. If we 

know G{ki,...,k]y) in Ri, by calculating the Fourier 
transform in Eq. we find the wave function ipB,c at 
some finite time t > in Ri. In Ref. [l^, -tpp^ki, . . . , /cat) 
was found to be proportional to the projection coeffi- 
cients . . . , fc^v) of the initial bosonic wave functions 
onto the LL eigenstates in free space. This connection 
allows us to calculate G(fci, . . . , kpf) for a few particles, 
which are in the ground state of a LL gas in a box po- 
tential and to study free expansion from such an 
initial state; this is performed in Sec. IIVI and Appendix 
[K\ For clarity, it should be noted that G(/ci, . . . , kjy) also 
depends on the coordinates Xj through the sgn(.Tj — Xi) 
terms {G(ki, . . . , fcjv) is not the Fourier transform of the 
bosonic wave function [III). 

In principle, from the time-dependent LL wave func- 
tion iIjb,c{xi, . ■ . ,XN,t) one can extract the physically 
relevant observables (in practice, this is a difhcult task). 
Here we consider one-body observables contained within 
the reduced single-particle density matrix (RSPDM), 



III. ASYMPTOTIC FORM OF THE 
MOMENTUM DISTRIBUTION 

In this section wc derive the asymptotic form of 
the momentum distribution of a Lieb-Liniger gas af- 
ter free expansion from an initially localized state de- 
fined by G(fci, . . . , fc^v) [we should keep in mind that 
G(fci, . . . , /cat) also depends upon the coordinates Xj via 
the sgn functions, see Eq. (jlOp j. The momentum distri- 
bution defined in Eq. (|12[) can be rewritten by using Eqs. 
(HD and mi) as 

nsik) ^ ^ j dxdye*''^^"^^ j dx2---dxN 

X f f dfci---(ifcjvG(fci,...,fcAr)e*^f=i('=^^^-'=?*)') 

X ffdq,- • •dgjvG(gi, . . . ,gjv)e'^" 

N f f 
~ — / dxi ■ ■ ■ dxN / dxdy dki ■ ■ ■ dkpj dqi ■ ■ ■ dq^ 
27r 7 J 

X G{ki,...M)*\x,=xG{qi,...,qN%,=ye'^, (14) 
where the phase </> is 

N 

(j){ki,. . . ,fcjv,9i, ■ ■ ■ ,qN,x,y) = - ^{kjXj - fc|t) 

N 

+ ^^{qjXj - g|t) - kix + kft + qiu - qft + k{x - y). 
i=2 







d4> 


d(j) 


d(j> 




d4> 


dxNtpB.ci^, X2, ■ 


■ ,XN,t)* 


dkj 


k'^ ~ 9qj 


1 dx 

1, 


X 


dy 



X 1pBAy^^2, ■ ■ ■,XN,t). 



(11) 



The SP density in real space is simply pB,cix,x,t), 
whereas the momentum distribution is defined as 

nB{k,t)^^ j dxdye'''^^-y^PB,c{x,y,t). (12) 

The eigenfunctions of the RSPDM, $i(a;, t) are called the 
natural orbitals (NOs), 



dxpBA^^V^t)'^i{^^i) A,(t)$,(?/,t), i = 1,2, . . . ; 

(13) 

the eigenvalues \i{t) are the occupancies of these orbitals. 
Apparently, in a noncquilibrium situation, the effective 
single particle states <^i{x,t) and their occupancies \i{t) 
may change in time. 



The integrals over ki, . . . ,kN,qi, . . . ,qN,x,y in Eq. (|T4l) 
are evaluated with the stationary phase approximation. 
The point of stationary phase is defined by the following 
equations: 



= 0, for 1 < j < A^. 



The stationary phase point is: 

k'j = q'j = x.j/{2t), for 2 < j < A^, 
ki = qi = k, and 
X ^ 2kt. 

The phase (j) can be rewritten as 

N 



(15) 



t 



''1 \2 



i=2 



2t 



2t' . 



4- [(fc -ki)x + kit] - [{k - qi)y + qlt]. 



(16) 



We notice that (j){k^, . . . ,kj^,q-y, . . . ,qj^,x ,y ) = 0. In the 
stationary phase approximation, the function G in Eq. 
is evaluated at the stationary phase point defined in 
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Eq. ([15]), which yields 
N 



nB,oo{k) 



2ti 



dx2 ■ ■ ■ dx 



N 



G{k. 



X2_ 

2f 



■' 2t 



N 



N-l 



dq2e 



(27r)3 



dx2 ■ ■ ■ dx 



N 



G{k. 



X2_ 

2i' 



N-1 



-it{q2- 



2< • 



JV-1 



(17) 



It is convenient now to introduce variables — Xj/t; 
from (|17p we obtain the asymptotic form of the momen- 
tum distribution of a freely expanding LL gas 

J d^2---d^N |G(fc,6/2,...,eiv/2)p. (18) 

We note that in the asymptotic regime, the momentum 
distribution acquires the same functional form as the 
asymptotic SP density. In the asymptotic regime, the SP 
density exhibits self-similar (ballistic) expansion (this is 
not true in the transient period preceding the asymptotic 
regime, sec Rcf. [l3|). It is most convenient to express 
the asymptotic SP density in variable £^ = x/t (see Ref. 

[F 



Poo(e)« J d^2---d^N |G(e/2,6/2,...,eiv/2)|2; (19) 

we normalize Poo(C) such that / Poo(C)^C = The vari- 
able = x/t has units of velocity; the self-similar asymp- 
totic SP density can be interpreted as the distribution 
of velocities of particles in a gas, which is in a simple 
manner related to the momentum distribution nB,oo{k). 

Equation can be thought of as a generalization of 
the dynamical fermionization of the momentum distribu- 
tion which has been demonstrated for a freely expanding 
TG gas (c oo) [l^, |T^. Free expansion in the TG 
regime is solved by the Fermi-Bose mapping [l|, [2^ . In 
this regime, the SP density is identical on both sides of 
the map. Since fermions are noninteracting, the asymp- 
totic form of the SP density (for both TG bosons and 
free fermions) is identical to the fermionic momentum 
distribution, which does not change in time. Equations 
(fT5)) and (fTO|) immediately yield that the asymptotic mo- 
mentum distribution for TG bosons has the same shape 
as the asymptotic SP density, which has the shape of 
the fermionic momentum distribution, i.e., we obtain the 
result of Refs. [l^,[lB|. We also note that equivalent rela- 
tion between the asymptotic SP density and momentum 
distribution was found in Ref. [g^l for a different model 
with emphasis that the time of flight measurements do 
not give the initial momentum distribution. The derived 
formula (jl8p is verified numerically on a particular ex- 
ample in the next section. 



IV. FREE EXPANSION FROM A BOX: 
DYNAMICS OF THE MOMENTUM 
DISTRIBUTION AND THE OCCUPANCIES \{t) 

In this section we calculate free expansion of three LL 
bosons, which are initially (at < = 0) in the ground state 
in an infinitely deep box of length L = tt. The analytical 
expression for the LL box ground state has been found 
in Ref. [l^. By using this result it is straightforward 
to calculate G(fci, k2, k^) (which depends on the interac- 
tion strength c) for this particular initial condition; we 
have outlined this calculation in Appendix |^ for N par- 
ticles. The next step is calculation of the Fourier integral 
in Eq. ^ , which is performed numerically by employing 
the Fast Fourier Transform algorithm. From the numer- 
ically obtained LL wave function 'ipB,c{xi,X2,X3,t) we 
calculate the momentum distribution nB{k,t), the SP 
density pc{x,t), natural orbitals and their occupancies, 
and study their evolution during free expansion from the 
box ground state. 

First let us explore the dynamics of the wave function 
''pB,cixi, X2, X3,t). Figure [T] displays contour plots of the 
probability density \4>B,c(fi,X2,X3,t)\^ for c = 1, at two 
different times, t = and t = 3. We see that as the 
LL gas expands, the probability density decreases at the 



hyperplanes Xi 



{i ^ j) where the particles are in 



contact. This is in agreement with the result of Ref. [13j . 
where it was shown (by using the stationary phase ap- 
proximation) that the leading term of '4>B,c{£,it, £,2t, ^st, t) 
has Tonks- Girardeau form for sufficiently large t; that is, 
the leading term is zero for = (i ^ j). However, this 
does not necessarily mean that the properties of such an 
asymptotic state correspond to the properties of a TG 
gas, which was usually studied in the ground state of 
some external potential. For example, suppose that the 
initial state is a weakly correlated ground state in the 
box; despite the fact that, during expansion, the parti- 
cles get strongly correlated in the close vicinity of the 
hyperplanes of contact, the absence of correlations in the 
initial state survives as an overall feature th roug h to the 
asymptotic state (see the discussion in Ref. [1^ and the 
second item of Ref. [l^l)- Thus, even though that the 
asymptotic state is described by a wave function with the 
TG structure, the physical properties of the expanded gas 
can considerably differ from the properties of a TG gas. 

In order to further study the properties of the state in 
expansion. Fig. [2] illustrates the occupation of the low- 
est natural orbital in time, Ai(<), for several values of 
c. The asymptotic values of the occupancies, which are 
obtained by using the asymptotic forms of the wave func- 
tions [l^, are indicated with horizontal lines. We observe 
that the occupancy of the leading NO, Ai(t), decreases 
during time evolution. However, for the plotted interac- 
tion strengths, the decrease of Xi{t) is not too large. This 
means that the coherence of the system (described by the 
occupations of the natural orbitals) for the plotted pa- 
rameters only partially decreases during free expansion 
due to the interactions. It should be noted that in the 




FIG. 2: (color online) The lowest natural orbital \i{t) as a 
function of time for three values of c. Red diamonds (dashed 
line) is for c = 0.25, black circles (solid line) for c — 1, and 
blue squares (dotted line) for c = 5; the lines connecting the 
markers are guides for the eye. The corresponding horizon- 
tal lines without markers denote the asymptotic occupancies, 
calculated from the asymptotic wave functions [l^ . 
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FIG. 1: (color online) Contour plots of |'!/'s,c(0, X2,X3, t)\'^ for 
c = 1 at (a) t=0, and (b) t=3. As the time t increases, the 
probability density at the hyperplanes where particles are in 
contact decreases. 



TG limit c — > cxd, for hard-core bosons on the lattice 
it has been shown that the leading natural orbitals 
slightly increase during free expansion [l5[ , which differs 
from the finite c results obtained here. It is reasonable 
to associate the decrease of Ai {t) to the change of the LL 
wave functions at the hyperplanes of contact; this change 
does not occur in the TG regime, where the wave func- 
tions are zero at the contact hyperplanes at any time of 
the expansion. 

Let us explore the dynamics of the momentum dis- 
tribution 77,5 (fc,t). and its connection to the SP density 
Pc{x,t) at large times t. The time-evolution of pc{x,t) 
and nB{k,t) is illustrated in Figs. [3] and SI we display 
X- and fc-space densities for various values of the param- 
eter c, at several times t. Initially, all momentum dis- 
tributions have a typical bosonic property: they peak at 
fc = 0. We observe that the qualitative changes in the 
shape of nB{k,t) are more pronounced for larger values 



of c. Circles in Figs. [3] and S] show the asymptotic values 
calculated by using Eqs. p8|) and (fT9|l . We see that at 
the maximal value of time t in the plots, the momentum 
distribution agrees well with that obtained with the sta- 
tionary phase approximation in Eq. (|18|) . Our numerical 
calculation is in agreement with the findings presented in 
Eqs. (fTSl) and We would like to point out that, even 
though the obscrvablcs nB{k,t) and pc{x,t) are well ap- 
proximated by the stationary phase approximation at the 
maximal expansion time reached in our numerical simu- 
lations (see Figs. [3]and|3]), the system is strictly speaking 
not yet fully in the asymptotic regime (e.g., note that the 
occupancies of the natural orbitals have not reached their 
asymptotic values) and even better agreement should be 
expected at larger times. Unfortunately, the maximal 
time allowed in our numerical calculations is limited by 
the computer memory and time. 

In order to further study the asymptotic forms of the 
momentum distribution and the SP density, let us calcu- 
late the asymptotic expansion velocity as a function of 
the interaction parameter c. Since different parts of the 
cloud expand at different velocities, a definition of this 
quantity has a certain degree of freedom. Here we define 
this quantity as a root mean square of the asymptotic SP 
density [i^ in variable ^ — x/t (i.e., velocity): 



^^ = ]J^ JePooiOdC, (20) 

the factor 1/A^ simply reflects the fact that Poo(0 is nor- 
malized to the number of particles N. The asymptotic 
velocity is connected to the total energy E stored 
in the system. During free expansion, the interaction 
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FIG. 3: (color online) Evolution of the 2;-space density in 
time for various interaction strengths c: (a) c = 0.25, at t = 
(red dotted line), t — 2 (solid black line), t = 4 (blue dashed 
line); (b) c = 1, at t = (red dotted line), t = 2 (solid black 
line), t = 4 (blue dashed line); (c) c — 10, a,t t — (red 
dotted line), t = 1 (solid black line), t = 3 (blue dashed line). 
The asymptotic a:-space density Poo(C) (circles), is plotted as 
a function of x = corresponding to the largest time in each 
subplot. 



FIG. 4: (color online) Evolution of the momentum distribu- 
tion in time for various interaction strengths c. The lines and 
colors for different values of c and t are identical as in Fig [3] 
Solid black and blue dashed line are almost indistinguishable. 



By using Eqs. ^TE\\ and p^ . we obtain 



Co 



(22) 



energy is transferred to the kinetic energy; in the asymp- 
totic regime all of the energy is kinetic, and it can be 
expressed via the momentum distribution: 



\k)dk. 



(21) 



that is, E = N^1^/A which is the classical expression for 
the kinetic energy of N particles with velocity ^oo and 
mass 1/2 (recall that we use units where the kinetic en- 
ergy operator in Eq. ^ is — X]i=i ^^/^^?)- The quan- 
tities \/~E and are displayed in Fig. O for various 
values of the interaction strength c; the plots underpin 
Eq. (|^^ . The total energy was calculated simply as 
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5 10 

c 

FIG. 5: (color online) Asymptotic expansion velocity, ^oo 
(squares, dashed line), and the square root of the total energy, 
(circles, solid line) for various interaction strengths c; 
lines serve to guide the eye (see text for details). 



the expansion velocity in asymptotic regime and pointed 
out its relation to the overall energy of the system. 

In order to gain further understanding of a freely ex- 
panding LL gas, it would be desirable to investigate tran- 
sient dynamics of the observables for larger number of 
particles, and also for different initial conditions (e.g., 
the ground state of a LL gas in different initial trapping 
potentials). 
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APPENDIX A: THE FUNCTION G(fci, . . . , fcjv) FOR 
THE BOX GROUND STATE 



E = qf + q2 + q'^ where quasimomenta qi are obtained 
by solving transcendental Bethe equations for the initial 
state (see Appendix!^. The asymptotic velocity was 
obtained via Eq. (|21[) by numerical integration. Our nu- 
merical calculations are in good agreement (better than 
99%) with Eq. (|22p : we attribute the discrepancy to in- 
accuracy of the numerical integration. 



V. CONCLUSION 

We have numerically studied free expansion of a few 
Lieb-Liniger bosons, which are initially in the ground 
state of an infinitely deep hard-wall trap. This numerical 
calculation has been carried out by employing a standard 
Fourier transform, as follows from the Fermi-Bose trans- 
formation for a time-dependent Lieb-Liniger gas. We 
have studied the evolution of the momentum distribution, 
the real-space single-particle density, and the occupan- 
cies of natural orbitals, both in the non-trivial transient 
regime of the expansion and asymptotically. We have 
derived analytically (by using the stationary phase ap- 
proximation) the formula which connects the asymptotic 
shape of the momentum distribution and the initial state. 
For sufficiently large times the momentum distribution 
coincides (up to a scaling transformation) with the shape 
of the real-space single-particle density (the expansion is 
asymptotically ballistic). This result can be considered 
as a generalization of the dynamical fermionization of the 
momentum distribution in the Tonks- Girardeau regime, 
which has been pointed to occur in the course of free ex- 
pansion [l^, . We have shown that the occupancy of 
the lowest natural orbital of the system decreases with 
time while approaching its asymptotic value. This was 
related to the build-up of correlations of the hyperplanes 
of contact of the particles. Finally, we have calculated 



In Sec. IIVI we have studied free expansion of three 
LL bosons, which are initially (at t = 0) in the ground 
state in an infinitely deep box of length L = n. Here we 
present exact analytical expression for function G{{k}) = 
G{ki, . . . , kpf) for this particular case. First, we use the 
connection between '0i^(fci, . . . , fc^v) and the projection 
coefficients 6(fci, . . . , A: a?) of the initial bosonic wave func- 
tions onto the LL eigenstates in free space [l^ to rewrite 
the expression for G: 



G{{k}) = NW{{k})b{{k}) 



X Yi [sgn(xj 

l<i<j<N 



,) + -(fc, -fc,)]. (Al) 



Here, J\f{{k}) is the normalization constant for LL eigen- 
states in free space jssj . 



(27r)^Ar![] 



(A2) 



N{{k}) \ 

and coefficients 6({fc}) are found by using the solution 
for the LL box ground state [2^ , 

h{{k}) (X AA({fc}) ^(-1)^' n [i - ^( V, - V.) 



l<i<j<N 

ic 



qi + <ij 



rL/2 i-L/2 r-Ljl 

X / dx\ / dxi • • • / 

-L/2 J x\ Jxj^^x 



)(1 



qpi — <iPj 



-) 




qpj - kp,,j\ Xj 



L 



qpj- 



(A3) 



In the expression above, summations are taken over all 
permutations P and P which are of order N, whereas 
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the set {e} is defined such that each is either +1 or — 1 
(here i ^ 1, . . . ,N, i.e., there are 2^ combinations in the 
set {e}). The ground state quasimomenta are defined as 
Qi = for I = 1,...,A^, and their magnitudes \qi\ 

are found by solving (numerically) the system of coupled 
transcendental equations [l^l 

= TT + V ( tan-i , , ^ , , + tan"^ , , ^ , , ) • 

(A4) 



Finally, let us mention that the constant of proportional- 
ity in Eq. (|A3p is fixed such that the wave function iJjb^c 
is properly normalized. 
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